home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DGELB.FOR next >
Text File  |  1985-11-29  |  7KB  |  249 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DGELB
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH A
  8. C           COEFFICIENT MATRIX OF BAND STRUCTURE.
  9. C
  10. C        USAGE
  11. C           CALL DGELB(R,A,M,N,MUD,MLD,EPS,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           R      - DOUBLE PRECISION M BY N RIGHT HAND SIDE MATRIX
  15. C                    (DESTROYED). ON RETURN R CONTAINS THE SOLUTION
  16. C                    OF THE EQUATIONS.
  17. C           A      - DOUBLE PRECISION M BY M COEFFICIENT MATRIX WITH
  18. C                    BAND STRUCTURE (DESTROYED).
  19. C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
  20. C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
  21. C           MUD    - THE NUMBER OF UPPER CODIAGONALS (THAT MEANS
  22. C                    CODIAGONALS ABOVE MAIN DIAGONAL).
  23. C           MLD    - THE NUMBER OF LOWER CODIAGONALS (THAT MEANS
  24. C                    CODIAGONALS BELOW MAIN DIAGONAL).
  25. C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS
  26. C                    RELATIVE TOLERANCE FOR TEST ON LOSS OF
  27. C                    SIGNIFICANCE.
  28. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  29. C                    IER=0  - NO ERROR,
  30. C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
  31. C                             TERS M,MUD,MLD OR BECAUSE OF PIVOT ELEMENT
  32. C                             AT ANY ELIMINATION STEP EQUAL TO 0,
  33. C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  34. C                             CANCE INDICATED AT ELIMINATION STEP K+1,
  35. C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
  36. C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
  37. C                             ABSOLUTELY GREATEST ELEMENT OF MATRIX A.
  38. C
  39. C        REMARKS
  40. C           BAND MATRIX A IS ASSUMED TO BE STORED ROWWISE IN THE FIRST
  41. C           ME SUCCESSIVE STORAGE LOCATIONS OF TOTALLY NEEDED MA
  42. C           STORAGE LOCATIONS, WHERE
  43. C             MA=M*MC-ML*(ML+1)/2    AND    ME=MA-MU*(MU+1)/2    WITH
  44. C             MC=MIN(M,1+MUD+MLD),  ML=MC-1-MLD,  MU=MC-1-MUD.
  45. C           RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
  46. C           IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN SOLUTION
  47. C           MATRIX R IS STORED COLUMNWISE TOO.
  48. C           INPUT PARAMETERS M, MUD, MLD SHOULD SATISFY THE FOLLOWING
  49. C           RESTRICTIONS     MUD NOT LESS THAN ZERO
  50. C                            MLD NOT LESS THAN ZERO
  51. C                            MUD+MLD NOT GREATER THAN 2*M-2.
  52. C           NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
  53. C           RESTRICTIONS ARE NOT SATISFIED.
  54. C           THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
  55. C           PARAMETERS ARE SATISFIED AND IF PIVOT ELEMENTS AT ALL
  56. C           ELIMINATION STEPS ARE DIFFERENT FROM 0. HOWEVER WARNING
  57. C           IER=K - IF GIVEN - INDICATES POSSIBLE LOSS OF SIGNIFICANCE.
  58. C           IN CASE OF A WELL SCALED MATRIX A AND APPROPRIATE TOLERANCE
  59. C           EPS, IER=K MAY BE INTERPRETED THAT MATRIX A HAS THE RANK K.
  60. C           NO WARNING IS GIVEN IF MATRIX A HAS NO LOWER CODIAGONAL.
  61. C
  62. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  63. C           NONE
  64. C
  65. C        METHOD
  66. C           SOLUTION IS DONE BY MEANS OF GAUSS ELIMINATION WITH
  67. C           COLUMN PIVOTING ONLY, IN ORDER TO PRESERVE BAND STRUCTURE
  68. C           IN REMAINING COEFFICIENT MATRICES.
  69. C
  70. C     ..................................................................
  71. C
  72.       SUBROUTINE DGELB(R,A,M,N,MUD,MLD,EPS,IER)
  73. C
  74. C
  75.       DIMENSION R(1),A(1)
  76.       DOUBLE PRECISION R,A,PIV,TB,TOL
  77. C
  78. C     TEST ON WRONG INPUT PARAMETERS
  79.       IF(MLD)47,1,1
  80.     1 IF(MUD)47,2,2
  81.     2 MC=1+MLD+MUD
  82.       IF(MC+1-M-M)3,3,47
  83. C
  84. C     PREPARE INTEGER PARAMETERS
  85. C        MC=NUMBER OF COLUMNS IN MATRIX A
  86. C        MU=NUMBER OF ZEROS TO BE INSERTED IN FIRST ROW OF MATRIX A
  87. C        ML=NUMBER OF MISSING ELEMENTS IN LAST ROW OF MATRIX A
  88. C        MR=INDEX OF LAST ROW IN MATRIX A WITH MC ELEMENTS
  89. C        MZ=TOTAL NUMBER OF ZEROS TO BE INSERTED IN MATRIX A
  90. C        MA=TOTAL NUMBER OF STORAGE LOCATIONS NECESSARY FOR MATRIX A
  91. C        NM=NUMBER OF ELEMENTS IN MATRIX R
  92.     3 IF(MC-M)5,5,4
  93.     4 MC=M
  94.     5 MU=MC-MUD-1
  95.       ML=MC-MLD-1
  96.       MR=M-ML
  97.       MZ=(MU*(MU+1))/2
  98.       MA=M*MC-(ML*(ML+1))/2
  99.       NM=N*M
  100. C
  101. C     MOVE ELEMENTS BACKWARD AND SEARCH FOR ABSOLUTELY GREATEST ELEMENT
  102. C     (NOT NECESSARY IN CASE OF A MATRIX WITHOUT LOWER CODIAGONALS)
  103.       IER=0
  104.       PIV=0.D0
  105.       IF(MLD)14,14,6
  106.     6 JJ=MA
  107.       J=MA-MZ
  108.       KST=J
  109.       DO 9 K=1,KST
  110.       TB=A(J)
  111.       A(JJ)=TB
  112.       TB=DABS(TB)
  113.       IF(TB-PIV)8,8,7
  114.     7 PIV=TB
  115.     8 J=J-1
  116.     9 JJ=JJ-1
  117. C
  118. C     INSERT ZEROS IN FIRST MU ROWS (NOT NECESSARY IN CASE MZ=0)
  119.       IF(MZ)14,14,10
  120.    10 JJ=1
  121.       J=1+MZ
  122.       IC=1+MUD
  123.       DO 13 I=1,MU
  124.       DO 12 K=1,MC
  125.       A(JJ)=0.D0
  126.       IF(K-IC)11,11,12
  127.    11 A(JJ)=A(J)
  128.       J=J+1
  129.    12 JJ=JJ+1
  130.    13 IC=IC+1
  131. C
  132. C     GENERATE TEST VALUE FOR SINGULARITY
  133.    14 TOL=EPS*PIV
  134. C
  135. C
  136. C     START DECOMPOSITION LOOP
  137.       KST=1
  138.       IDST=MC
  139.       IC=MC-1
  140.       DO 38 K=1,M
  141.       IF(K-MR-1)16,16,15
  142.    15 IDST=IDST-1
  143.    16 ID=IDST
  144.       ILR=K+MLD
  145.       IF(ILR-M)18,18,17
  146.    17 ILR=M
  147.    18 II=KST
  148. C
  149. C     PIVOT SEARCH IN FIRST COLUMN (ROW INDEXES FROM I=K UP TO I=ILR)
  150.       PIV=0.D0
  151.       DO 22 I=K,ILR
  152.       TB=DABS(A(II))
  153.       IF(TB-PIV)20,20,19
  154.    19 PIV=TB
  155.       J=I
  156.       JJ=II
  157.    20 IF(I-MR)22,22,21
  158.    21 ID=ID-1
  159.    22 II=II+ID
  160. C
  161. C     TEST ON SINGULARITY
  162.       IF(PIV)47,47,23
  163.    23 IF(IER)26,24,26
  164.    24 IF(PIV-TOL)25,25,26
  165.    25 IER=K-1
  166.    26 PIV=1.D0/A(JJ)
  167. C
  168. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
  169.       ID=J-K
  170.       DO 27 I=K,NM,M
  171.       II=I+ID
  172.       TB=PIV*R(II)
  173.       R(II)=R(I)
  174.    27 R(I)=TB
  175. C
  176. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN COEFFICIENT MATRIX A
  177.       II=KST
  178.       J=JJ+IC
  179.       DO 28 I=JJ,J
  180.       TB=PIV*A(I)
  181.       A(I)=A(II)
  182.       A(II)=TB
  183.    28 II=II+1
  184. C
  185. C     ELEMENT REDUCTION
  186.       IF(K-ILR)29,34,34
  187.    29 ID=KST
  188.       II=K+1
  189.       MU=KST+1
  190.       MZ=KST+IC
  191.       DO 33 I=II,ILR
  192. C
  193. C     IN MATRIX A
  194.       ID=ID+MC
  195.       JJ=I-MR-1
  196.       IF(JJ)31,31,30
  197.    30 ID=ID-JJ
  198.    31 PIV=-A(ID)
  199.       J=ID+1
  200.       DO 32 JJ=MU,MZ
  201.       A(J-1)=A(J)+PIV*A(JJ)
  202.    32 J=J+1
  203.       A(J-1)=0.D0
  204. C
  205. C     IN MATRIX R
  206.       J=K
  207.       DO 33 JJ=I,NM,M
  208.       R(JJ)=R(JJ)+PIV*R(J)
  209.    33 J=J+M
  210.    34 KST=KST+MC
  211.       IF(ILR-MR)36,35,35
  212.    35 IC=IC-1
  213.    36 ID=K-MR
  214.       IF(ID)38,38,37
  215.    37 KST=KST-ID
  216.    38 CONTINUE
  217. C     END OF DECOMPOSITION LOOP
  218. C
  219. C
  220. C     BACK SUBSTITUTION
  221.       IF(MC-1)46,46,39
  222.    39 IC=2
  223.       KST=MA+ML-MC+2
  224.       II=M
  225.       DO 45 I=2,M
  226.       KST=KST-MC
  227.       II=II-1
  228.       J=II-MR
  229.       IF(J)41,41,40
  230.    40 KST=KST+J
  231.    41 DO 43 J=II,NM,M
  232.       TB=R(J)
  233.       MZ=KST+IC-2
  234.       ID=J
  235.       DO 42 JJ=KST,MZ
  236.       ID=ID+1
  237.    42 TB=TB-A(JJ)*R(ID)
  238.    43 R(J)=TB
  239.       IF(IC-MC)44,45,45
  240.    44 IC=IC+1
  241.    45 CONTINUE
  242.    46 RETURN
  243. C
  244. C
  245. C     ERROR RETURN
  246.    47 IER=-1
  247.       RETURN
  248.       END
  249.